Robust quantum searching with spontaneously decaying qubits 



O 

o 

(N 

(N 

> 

^ ■ 

o ■ 

On : 

O . 

\0 ■ 

o ■ 

:^ 

Oh! 

> ■ 

■ 

a • 

^ : 



Robert J. C. Spreeu'vJB and Tom W. Hijmans 

Van der Waals-Zeeman Institute, University of Amsterdam, 
Valckenierstraat 65, 1018 XE Amsterdam, The Netherlands 
(Dated: February 1, 2008) 

We present a modification of the standard singie-item quantum searcfi procedure tfiat acquires 
robustness from spontaneous decay of the qubits. This damps the usual oscillation of populations, 
driving the system to a steady state with a strongly enhanced population of the solution. Numerical 
evaluation of the steady state was performed for up to 36 qubits. The huge size of the state space in 
our analysis is dealt with by exploiting a symmetry in the master equation that reduces the scaling 
of computer resources from exponential to polynomial. Based on these results we estimate that 
an error-free solution can be retrieved from the steady state after 0(log log A'') repetitions, with 
near-unit probability. This brings the overall scaling to 0{^/N log log A'^), only slightly worse than 
for the ideal quantum case. 

PACS numbers: 03.67.Pp, 03.67.Lx 



I. INTRODUCTION 

The major challenge in quantum information science 
and technology is to achieve complete control over quan- 
tum systems. In this light it has been realized early on 
that loss of coherence can be a prohibitive problem for 
quantum computers P, [3j 111- In general, spontaneous 
and incoherent processes are therefore best avoided where 
ossible, or counteracted by quantum error correction 
Q. From this perspective it is remarkable that in 
specific cases dissipation can also be put to use. Earlier 
work has shown that dissipation and noise can assist in 
the production of, or even generate entanglement [6, 7, 8]. 
Here we modify Grover's quantum search algorithm 9] to 
yield a different, more robust, operating procedure in the 
presence of local decay. Instead of applying the search it- 
erations a prescribed number of times, we let the system 
relax to a steady state by spontaneous decay. 

Several authors have investigated the effects of errors 
on the performance of quantum algorithms. Most stud- 
ies have concentrated on unitary errors, including ran- 
dom errors (noise) or systematic (static) imperfections 
[lOl [ill , m, [H, In this paper we study errors 

of a nonunitary, dissipative nature. These have pre- 
viously been investigated in various different contexts 
P, 0, H, [3 • Experimentally, the effects of losses were 
studied in an optical, classical-wave analog of quantum 
searching [TgI. IT7|. 

Robustness against errors and resistance to decoher- 
ence has previously been reported for adiabatic quantum 
computation p^ . Although aiming for similar benefits, 
a few essential differences with our approach presented 
here are worth pointing out. Adiabatic quantum com- 
putation requires dynamic control over the Hamiltonian 
as the computation progresses. In our approach the dy- 
namic control is entirely absent, to the extent that the 



'Electronic address: ' spreeuwQscience . uva. nl] 



initialization as well as the timing of the computation are 
rendered superfluous. A second essential difference con- 
cerns the role of the dissipation. Whereas adiabatic quan- 
tum computation has been reported to be robust against 
decoherence, in our modified quantum search procedure a 
specific type of dissipation plays an active role: it is used 
to drive the system to a steady state. Consequently, the 
acquired robustness is of a different nature. 

The idea of the present paper can be understood as fol- 
lows. The quantum search algorithm essentially drives a 
rotation in a two-dimensional subspace of the full Hilbert 
space, spanned by the initial state and the solution. Noise 
and incoherent processes will usually lead to loss of per- 
formance by "leakage" of population out of this subspace 
[l3| . We show that we can prevent the system from stray- 
ing too far from this two-dimensional subspace, if the 
noise is of a "natural" origin: spontaneous qubit decay. 

We numerically simulate a single-item quantum search 
for a marked solution, in the presence of spontaneous de- 
cay, on our classical computer. In principle we would 
suffer from an exponential scaling of resources [l^. For- 
tunately, for the problem at hand we identify a symmetry 
that reduces the scaling to polynomial in the number of 
qubits, with ^ density matrix elements. By exploit- 
ing this symmetry we can handle up to 36 qubits on a 
desktop personal computer. 

The remainder of the paper is structured as follows. In 
Sec. In] we adopt the Hamiltonian description of a single- 
item quantum search and we extend this description to 
include also decay processes. We give the resulting mas- 
ter equation for the density matrix that includes our spe- 
cific choice for the decay model, which is crucial for the 
results in this paper. In Sec. IIIII we identify a symmetry 
in the steady state that allows us to drastically reduce the 
amount of data in a reduced density matrix. In Sec. IIVI 
we discuss the numerical results obtained by solving for 
the steady state. We discuss how one would obtain an 
error-free solution by repetition and majority voting. 
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II. FORMAL DESCRIPTION 
A. Search Hamiltonian 

In the usual formulation of Grover's quantum search 
algorithm [9] , one first initializes a quantum register with 
all qubits in the state |0), and applies a bitwise Hadamard 
operation H^'^. Throughout this paper we indicate the 
number of qubits by q. The bitwise Hadamard yields 
a symmetric superposition of all basis states, (-ff|0))®*, 
which we write as 
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To this state one then repeatedly applies a unitary search 
iterator = (2|s)(s| — — 2\w){w\), where w is the 
solution, and / the identity operation. As a result, the 
population of the solution pww oscillates between and 
1. The readout is usually performed after (7r/4)ViV — 
(7r/4)2''/^ iterations, when p^j^ « 1. 

In this paper we shall use an alternative formulation 
in terms of time-continuous evolution ip = with 
the search Hamiltonian given by 



\w){w\ + \s){s\. 



(2) 



This formulation, which has been introduced by Farhi 
and Gutmann [20|], yields very similar behavior, with an 
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FIG. 1: (Color online) (a) Evolution of the population of 
the solution p^nm without decay (red dashed line) and with 
r = 0.03 (blue solid line), for q = 6 qubits. (b) Steady state 
populations p^x for all basis states for the above CEise with 
decay. The solution is a: = 0; equal colors indicate equal 
Hamming distance to the solution. 



oscillation period of tt^/N = tt 2'/^ . The Hamiltonian 
formulation has previousl y b een applied to discuss adia- 
batic quantum searching |2ll . [2^ . [23j . We choose it here 
because a time-continuous description is more convenient 
when determining steady states. Our conclusions remain 
valid also for the usual iterative version. We have verified 
this numerically for small qubit numbers q < 10, where 
iterations with the unitary were alternated with finite 
time intervals of qubit decay. For larger qubit numbers 
we expect the iterative algorithm to resemble a time- 
continuous process ever more closely. 



B. Description of the decay 

The oscillatory search behavior mentioned above is 
reminiscent of the Rabi oscillation in a two-level sys- 
tem coupled to a resonant driving field. In the latter 
case spontaneous decay damps the Rabi oscillation and 
relaxes the driven two-level system to a steady state. 
Guided by this analogy we now introduce spontaneous 
decay for the qubits, driving each qubit toward the state 
H\0) = (1/V2)(|0) + |1)) and the syste m as a whole to- 
ward the uniform superposition \s). 

This choice of |s) as the "ground state" of the decay 
is essential. In a loss- free search, starting with \s) as the 
initial state, the search Hamiltonian !Km) drives a rota- 
tion in a two-dimensional subspace spanned by \w) and 
|s). The quantum state is then confined to this subspace. 
Dissipation or decoherence will in general lead to leakage 
of population out of this subspace. In our case, however, 
the decay prevents the quantum state from straying too 
far from |s) and thus from the subspace. 

By our use of the word "decay" we mean to indicate 
that the process is unidirectional. We do not mean that 
|s), the state the system decays into, is the lowest eigen- 
state of In fact, all eigenvalues are zero, except the 
two highest ones, 1 ± l/ViV, corresponding to the eigen- 
states |s) ± |w) (not normalized). The choice to have the 
subspace {|s),|w)} at the highest energies rather than 
the lowest is arbitrary and unimportant for our results. 

It may seem unnatural to consider a process where 
the qubits decay to H\0) = (1/V2)(|0) + |1)), the natu- 
ral choice being perhaps decay to |0) (or |1)). However, 
our choice is equivalent to using the "natural" decay to 
the ground state |0), if the entire quantum search is per- 
formed in a basis transformed by For convenience, 
in our simulation we choose to let the qubits decay to 
H\0) and to perform the quantum search in the standard 
computational basis. 

The decay is described by means of a master equation 
for the density matrix p. We use a common model for 
spontaneous emission into a reservoir of modes at zero 
temperature. Thus using the Born-Markov approxima- 
tion, and including also the coherent Hamiltonian evolu- 
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tion, it has the following form (23 |. 

P^i[p, 'Kw] + ^'^{^Crpc\~ c\cr p - p C,^Q^ (3) 

4=0 

where F is a local (single-qubit) decay rate. The Ci are 
transition operators that describe the coupling to a reser- 
voir. They are chosen to be local lowering operators in 
the Hadamard-transformed basis, 

Ci^h®---® Hd~H (g) • • • (g) /2 (4) 

where the number of factors is q, and the single-bit low- 
ering operator d~ — |0)(1| appears at position i. At all 
other positions we have the 2x2 identity matrix 12- 

With the Ci thus defined, each qubit decays to the state 
H\0) = {1/V2){\0) + |1)). The collective "ground state" 
of the decay process is therefore the uniform superposi- 
tion \s). This is a crucial difference with previous work by 
Zhirov and Shepelyansky [l^ , where the decay drove the 
system to the state lO)®"^. In the latter case the buildup 
of population in \w) is not present in the steady state. 
In this paper we describe how this buildup does occur if 
the lowering operators Ci are defined in a basis that is 
rotated with respect to the standard computational ba- 
sis in which the coherent search process is defined. As 
mentioned above, we could have chosen not to rotate the 
dissipation basis and to rotate the search basis instead. 

C. Integration of the master equation 

We can now simulate the quantum search in the pres- 
ence of decay by integrating Eq. ([3]). An example is 
shown in Fig. [ija) for q — 6 qubits. For comparison 
we also show the result without decay. We notice that 
the oscillation of the population of the solution, pww, is 
damped by the decay. The steady state value is clearly 
much higher than 2~' « 0.016, which would be the value 
if the population were randomized entirely. A readout 
will thus reveal the solution with increased probability. 
The readout will typically contain errors in some fraction 
of bits. The full solution can be retrieved after repeating 
this procedure a modest number of times, as we show in 
SecHVl 

III. REDUCTION OF SCALING BY 
SYMMETRY 

We now investigate how the steady state population 
Pww depends on the loss rate F and how it scales with 
the number of qubits q. The latter is particularly nontriv- 
ial, since simulating a quantum computer on a classical 
one is usually inefficient. The number of entries in the 
density matrix is 2^"^ and quickly exhausts the resources 
of any classical computer. Fortunately we can exploit a 
symmetry in the problem to vastly reduce this number, 
and reduce the scaling to polynomial in q. 



This symmetry can be recognized in Fig. [Ijb) , show- 
ing how the population p^^ of any given basis state \x) 
depends only on the Hamming distance between x and 
w, defined as the number of bits in which x and w dif- 
fer. Inspection of the full density matrix reveals a similar 
symmetry for the off-diagonal elements, with p^y only a 
function of the Hamming distances d^wi dyw, and dxy 

To understand why this is the case we first perform 
a basis transformation that for a given solution w 
relabels the basis states as Su]\x) — \x') = \XOR{x,w)), 
where the XOR operation is to be taken bitwise. Note 
that Swlw) = |0) so that the Hamming distances become 
dxw — *■ dxo and similarly dyw — *■ dyo- The Hamming 
distance dxy is invariant under Sw 

The "shift" operator Sw transforms the search Hamil- 
tonian into that for w = 0: 

S^:K^SI^%o^\0){Q\ + \s){s\, (5) 

while it leaves the decay terms in the master equation 
invariant. This can be seen by writing = I2 ^ ■ ■ ■ ^ 
X (g) ■ ■ ■ ®) I2, where the number of factors is q and the 
Pauli matrix AT = |0)(1|-|-|1)(0| appears in every position 
i where the corresponding bit value of w is equal to 1. 
Defining c — Hd~H and c = XcX, one easily verifies 
that c^c = c^c and cp2C' = cp2C^ , for any arbitrary 2x2 
matrix p2- From this it follows that the decay terms in 
the master equation are unchanged by Sw 

From here on we consider the search problem with 
w = without loss of generality. This problem is sym- 
metric under bit swaps, i.e., under all possible permuta- 
tions of the qubits. This is immediately clear by inspec- 
tion: both |0) and |s) appearing in the search Hamilto- 
nian IKq are product states with all qubits in the same 
state. The decay terms in the master equation consist of 
sums of identical terms for all qubits and are therefore 
also invariant under qubit permutations. 

When searching for the steady state of the master 
equation, we can now restrict ourselves to density ma- 
trices that have bit swap symmetry. For a given ma- 
trix element pxy, it is convenient to introduce the bit 
pair counts noo, noi, «iOj a-nd nn, such that tiqo counts 
the number of instances where corresponding bits in x 
and y both have the value 0. The other bit-pair counts 
noi, niQ, and nn are defined similarly. Their relation to 
{dxo, dyQ, dxy) is 

dxo = nio + nil (6) 

dyo = uqi + nil (7) 

dxy = noi + niQ (8) 

q = nQo + noi + nio + nu (9) 

If we perform a permutation of qubits, so that x —i- x' 
and y ^ y' ^ corresponding bits in x and y are permuted 
in the same way. Therefore any permutation of qubits 
leaves the bit pair counts unchanged, so pxy = Px'y'- 

Instead of labeling the density matrix elements by in- 
dex pairs, we now label them by (noo, ^oii 't^iOj '^11): 

Pxy — Cnoo,noi ,nio,rill ■ (10) 
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FIG. 2: (Color online) Steady state population in the so- 
lution, Poo, as a function of the decay rate F. The data 
sets are for different numbers of qubits, from the top down: 
g = 6, 8, 10, 12, 16, 20, 24, 28. The decay rate has been scaled 
by the frequency of the quantum search oscillation, ~ 2"'''^. 

This is a crucial result: whereas x, y run from to 2"? — 1, 
the bit pair counts (or Hamming distances) run only from 
to q. Since the bit pair counts sum to q, the number of 
different combinations of bit pair counts is given by 

('^J^) = J('Z + l)(9 + 2)(g + 3). (11) 

This is a very modest number compared to the number 
of index pairs {x,y), equal to 2^*^. The mmibcr of dis- 
tinct entries in the density matrix is thus reduced tremen- 
dously. More importantly, the exponential scaling with 
q has been replaced by a polynomial one. For example, 
for q = 36, we have 2''^ k, 4.7 x 10^^ whereas bit swap 
symmetry reduces this to only 9139 distinct matrix ele- 
ments. 

We now reexpress the master Eq. (j3|) as the time evo- 
lution of the new density matrix, a. The expression is 
somewhat lengthy but straightforward and is given in 
the appendix. Using this equation we can now either 
integrate it in time, or directly solve a — Q to obtain 
the steady state. It should be noted however that we 
have hereby given up the possibility to start the time in- 
tegration with an arbitrary initial density matrix. The 
equation for a implicitly assumes the symmetry described 
above. The usual initial condition, starting from the su- 
perposition \s) does have this symmetry and can there- 
fore be simulated. 

If we solve directly for the steady state, we are not sub- 
ject to this limitation. The decay process ensures that the 
density matrix acquires the appropriate symmetry over 
time. The usual initialization step is thus unnecessary. It 
is worth noting that the time needed to develop bit swap 
symmetry is the same as the time to reach the steady 
state. This is not surprising, considering that both are 
driven by the same decay process. Numerically, this was 
observed by choosing a randomly chosen, (pure,) initial 
state /o, lacking bit swap symmetry. For obvious reasons 
this test is restricted to small numbers of qubits. 



From here on we solve only the equation (j = 0, which 
is a system of linear equations with as many variables. 
More precisely, these linear equations are not indepen- 
dent. This is resolved by replacing one of the linear 
equations by the unit trace condition Trp — 1, which 
translates to 

dM=0 ^ ^' ^ 

IV. NUMERICAL RESULTS FOR THE STEADY 
STATE 

In Fig. [2] we show the resulting steady state popula- 
tion in the solution, poo — fq,o,o,o, for various numbers 
of qubits q, and varying the loss rate F. The loss rate 
has been scaled by the frequency of the loss-free quan- 
tum search oscillation 2^''/^. We clearly recognize two 
different regimes, 29/2 r ^ ^ 29/2 F < 1. We find 
that in both regimes the steady state populations are 
independent of the decay rate F. There is a transition 
region around 2''/'^ F « 1 where the steady state popula- 
tion goes over from the low-F value to the high-F value. 
Not surprisingly, in the limit of strong decay the steady 
state approaches |s), so that all populations in the com- 
putational basis, pxx, approach the value 2"'. 

The low-F regime is thus reached when the single-qubit 
decay rate is well below the frequency of the loss-free 
search oscillation (the energy separation between the two 
highest eigenvalues of IH„). This means that the prob- 
ability that one particular qubit has decayed during one 
search period is small. However, the probability that at 
least one qubit decays during such a period can be high. 

On the basis of Fig. [2] we can now estimate the total 
search time. We see from this figure that for < 0.3 

we are essentially in the low F limit. The time T needed 
to reach the steady state is a few times F~^. For the 
example of Fig. [Ha) we have 2'^/'^T = 0.24 and FT = 3.6 
at T = 120. Thus the total search time is T « 15 x 2^/^, 
roughly ten times longer than the value of (7r/2)29/2 for 
the loss-free case. The above estimate states essentially 
that (i) the search oscillation must be well underdamped, 
and (ii) we need to wait for a few damping times. There- 
fore we expect the mentioned factor of ten to be essen- 
tially constant, i.e. not to scale with q. 

We now investigate the weak decay limit in some more 
detail. We choose, somewhat arbitrarily, 2''/^r = 0.005 
as representative of this low-F limit. The steady state 
population of the solution, poo, is then much larger than 
the other populations. We also see that poo decreases 
with increasing q. A plot of poo vs. 1/q is remarkably 
well fitted by a straight line. However this relationship 
is purely heuristical and extrapolates to negative popu- 
lation for q > 120. 

We define 29poo as an "amplification factor" , i.e. the 
factor by which poo is increased compared to the fully 
random case. This is also the ratio of the poo values 
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FIG. 3: Amplification factor 2' poo, in the limit of weak decay, 
versus the number of qubits q. 
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FIG. 4: (Color online) Distribution of steady state population 
over the Hamming distance d to the solution, in the weak de- 
cay limit. The different data sets are for different numbers 
of qubits, 5 = 6, 12, 18, 24, 30, 36. The insert shows the aver- 
age Hamming distance as a fraction of the number of qubits, 
versus q. 



in the low-F and the high-F limits. In Fig. [3] we show 
how the amplification varies with q in an approximately 
exponential fashion. 

It is interesting to compare the weak decay case to a 
strongly driven two-level system with spontaneous emis- 
sion. The populations of the ground and excited states 
then both approach the value 1/2. In the quantum search 
the population poo is generally much smaller, except for 
very small numbers of qubits. However, if we project 
the numerically obtained steady state on the subspace 
spanned by {|0), |s)}, we do find equal populations in 
these two states. 

The observation that poo becomes ever smaller for 
larger numbers of qubits may seem disappointing at first. 
However this is unjustly so. We now investigate how the 
remaining population 1 — poo is distributed and show that 
it is concentrated in states at small Hamming distance 
from the solution. 

The diagonal elements of the density matrix are re- 
trieved as pxx = cFq^d^oMfi,d^o- Sorting these by their 
Hamming distance d^o, and multiplying by their multi- 



plicity (a binomial coefficient) , we obtain the distribution 
of population over the Hamming distance d to the solu- 
tion. This is shown in Fig. m For small numbers of qubits 
the distribution is peaked at d = 0. For larger q the peak 
shifts away from zero. This can be understood qualita- 
tively because the multiplicity increases rapidly with d. 

For strong decay the mean of the distribution (d) = 
q/2. A readout then yields a wrong bit value for half 
the bits, on average. For weak decay we find that the 
"bit error rate" {d)/q is again almost constant, with a 
value of ^ = {d)/q « 0.28, as can be seen in Fig. ID The 
solution can now be obtained by repeating the search a 
limited number of times using a majority vote rule to 
decide about each individual bit value. We can estimate 
how many times R we must repeat the search and readout 
in order to increase the probability to find an error-free 
solution to above a preset probability. 

The probability that the majority vote yields the 
wrong result for a particular bit is obtained by sum- 
ming the binomial distribution, ( ^ ) ~ "i f^or 



n > {R + l)/2 (for odd R). For large R we can ap- 
proximate the summation using standard techniques and 
obtain an upper bound for the single-bit error rate after 
majority voting, 



1 (2v/e^ 



1 - 2C 



(13) 



The probability of an error in at least one bit out of q 
is then given by 1 — (1 — ^r)'^ < q£.R- For q = 29 and 
^ = 0.28 we find that q^R < 0.05 (0.01) for i? > 41 (53). 

It is important to note the scaling behavior oi^R. From 
Eq. (fT3|) we see that £,r decreases faster than exponential 
with R. Thus, if we demand that q^R drops below some 
preset error probability e, the number of repetitions will 
increase only as i? = 0(log(g/e)). The total search time 
then scales as 0(2fil'^ log((7/e)), somewhat longer than for 
the commonly studied relaxation-free case (2'^/^) but still 
much shorter than for a classical search (2'^). 



V. CONCLUSION 

In conclusion, we have shown that a modified, robust 
version of Grover's quantum search algorithm can be 
used if the qubits are subject to spontaneous decay. The 
search is conducted in a Hadamard transformed basis 
and the single-qubit decay rate must be smaller than the 
natural oscillation frequency of the search algorithm. 

A symmetry allowed us to numerically analyze the 
problem for up to 36 qubits on a standard desktop com- 
puter. This clearly opens up the possibility to investigate 
this problem for even larger number of qubits, at only 
polynomial costs. It remains an open question what will 
happen to the observed constant bit error rate of 0.28 
and the trend poo ~ 1/? in the limit q ^ oo. The lat- 
ter must break down at around 100 qubits. Possibly, an 
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analytic approximation in this limit could shed light on where 5i,j is the Kronecker delta and 31 and C are row 
this issue. and column sums, given by 

An obvious extension of the present work is multiple- 
item searching for which we expect similarly robust be- 
havior in the presence of relaxation. We can only spec- 
ulate on the question whether similar modifications may _^ ^ 
be applied to other quantum algorithms. ^ _ y-? f d^o \ ( ^ ~ 

joo=Ojii=0 ^ ^ ^ 

Acknowledgments '^joo,9— dxo— ^oojC^xo— *ii 

(A.2) 
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APPENDIX: RE-EXPRESSION OF THE 
MASTER EQUATION 

Using the symmetry described above, we relabel the 
density matrix elements Pxy using the bit pair counts, 
Pxy = ^naa,noi,nia.nii- The master cquation can then be 
re-expressed as follows, 

'-''"oo,"oi,nio,nii ~ 

2-«i(3^- e)+ 

r 

^ [ '^OO (f noo-l,noi,"io,"ii + l (^noo,noi,nio,nii) ^ 

nil (fnoo + l,"oi,"io,nii-l ~ f noo,"oi,"io,nii ) ^" 

noi {'^ <^noo,noi-l,nio,nii + l + 2 (7noo-|-l,noi -l,nio,raii " 

""^nooiraoi — l)"io-|-l,"ii ~ "'^nooi'ioi ,"io,"ii ) "I" 
f^lO (2 Crnoo,"oi,mo — l, "11-1-1 ''r ^noo + '^,noi,nio — l,nii~ 

~''^noo,"oi-|-l,"io — l,"ii ~ •'^"ooinoii'iioiraii ) ] 

(A.l) 
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